Delaunay Triangulations in Linear Time? 



We present a new and simple randomized algorithm for construct- 
ing the Delaunay triangulation using nearest-neighbor graphs for point 
location. Under suitable assumptions, it runs in linear expected time 
for points in the plane with polynomially bounded spread, i.e., if the 
ratio between the largest and smallest pointwise distance is polynomi- 
ally bounded. This also holds for point sets with bounded spread in 
higher dimensions as long as the expected complexity of the Delaunay 
triangulation of a sample of the points is linear in the sample size. 

Chan and Patra§cu [HI EJ presented o(iV log N) randomized algorithms for 
constructing Voronoi Diagrams of points in the plane (from which the De- 
launay triangulation can be computed in linear time and vice-versa) under 
suitable models of computation. Here we present an O(N) randomized algo- 
rithm for the Delaunay triangulation in the plane in a different model. The 
algorithm is not restricted to two dimensions and it runs in linear expected 
time as long as the expected complexity of the Delaunay triangulation of a 
random sample of the input points is linear in the sample size. An example 
of linear complexity Delaunay triangulation are suitably sampled (d — 1)- 
dimensional polyhedra in IRA Our algorithm locates points by combining 
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Algorithm 1: Incremental Construction using Nearest-Neighbor 
Graph 

Input: Finite point set P in IR d 
Output: Delaunay triangulation of P 

1 Split P into rounds Ri, . . . R m of doubling size with R\ of constant 

size, let Sj := Ui<i<j R i U = 1 • • ■ m )- 

2 Insert points in Ri into the Delaunay triangulation using history for 
point location. 

3 For k — 2, . . . , m insert points in R k into the Delaunay triangulation 
as follows: 

3.1 Set T fc _i «- R k , Ti «- (0 < i < k - 1), and j «- fc - 1. 

3.2 While Tj ^ and j > 0: compute NNG(I} U Sj) and from each 
connected component with no vertex in Sj add the first point of the 
component to then set j <— j — 1. 

3.3 While j < k — 1: locate 1} (if not empty) in DT(Sj+i) using history 
starting at DT(Sj); then locate T j+1 in DT(S , J+ i) by walking 
through the connected components starting at an already located 
point; then set j <— j + 1 

the history (i.e., the Delaunay tree [2113]) of a randomized incremental con- 
struction with a sequence of nearest-neighbor graph computations. For the 
nearest-neighbor graphs we use a recent result by Chan [7] that links well- 
separated pair decompositions to sorting. By the use of radix sort this results 
in a linear time algorithm for well-separated pair decompositions and as a 
consequence for nearest-neighbor graphs. We will use the same assumptions 
as Chan on the model of computation and the point set. The model of com- 
putation is a real-RAM with a constant-time restricted floor function that 
can be applied only if the resulting integer has 0(log N) bits. Restricting the 
floor function avoids issues about creating an unreasonably powerful model of 
computation. The input point set should have polynomially bounded spread, 
i.e., the ratio of the largest and smallest point to point distance should be 
bounded by a polynomial in the size of the point set. But also other combina- 
tions of models of computation and sorting algorithms can be used, resulting 
in a optimal running time asymptotically bounded by the time needed for 
sorting (see p] for details). 
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General Setup. We construct the Delaunay triangulation of a finite point 
set P C H d by a randomized incremental construction using a history. The 
point location is accelerated by locating points at intermediate levels in the 
history instead of the top, see Algorithm [TJ Given an insertion order we 
group the points into rounds R±, . . . ,R m in accordance with the order, i.e., 
the points in R t are in the insertion order before the points in Ri + \ for 
1 < i < m. The rounds double in size, i.e., is constant, and \R%+\ \ = 2\Ri\ 
(with possibly the exception of the last round for which \R m \ < 2i? m _i). Let 
Sj := Ui<i<j-^4 denote the points inserted in or before round j. Together 
with the history graph we store the Delaunay triangulations of the Sj. Note 
that the rounds are only used for facilitating the point location; the insertion 
order remains the same. 

Point Location in a Round. The points of the first round are located in 
the standard way using the history. At the beginning of round k (2 < k < m) 
the points of the round Rk are located in the Delaunay triangulation of S k -i 
using a family of sets Ti, . . . , T^-i (in every round a different family, thus 
more formally the family could be written as Tf, . . . , T^_ x ): Let T k _i := 
Rk- We compute the nearest-neighbor graph of T^-\ U S k -i- For connected 
components of the nearest-neighbor graph without a vertex in S k ~i we include 
the first point (according to the insertion order) of the component in a set 
Tk-2- We repeat the same procedure higher up in the history, i.e., we compute 
the nearest-neighbor graph of T fc __ 2 U for each connected component 

without a vertex in S k _ 2 we include the first point in a set T fe _ 3 , and so on. 
We stop this process with the construction of T (or earlier with Tj if is 
empty. For simplicity we describe the algorithm for the case that T is not 
empty). This yields a hierarchy of sets T fc _i D T k ^ 2 D ■ • O T . 

Now we locate the points in T in DT(S'i) by using the history, i.e., we 
use the history to find a conflicting simplex and then locally search for the 
simplex containing T . For locating T± we have the following situation: each 
connected component of the nearest-neighbor graph of T\ U S\ either has a 
vertex in Si or has a vertex in T , thus each component has a vertex already 
located in DT(S'i). We traverse each component starting at an already lo- 
cated vertex, e.g., by a depth first search. During the traversal we locate 
the traversed points in DT(S'i) by walking from an already located neighbor, 
i.e., we locally traverse the triangulation along the line segment between the 
two points. After locating the points in T\ in DT(Si) we locate them in 
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DT(S , 2 ) by using the history, starting not at the top of the history but at the 
simplices of DT(S'i) containing the points. By the same procedure we locate 
the points in T 2 in DT(S 2 ) and DT(S 3 ), and so on, until we have located the 
points in T k _i = R k in DT(Sfc_i). Finally we insert the points of R k into 
the Delaunay triangulation, where a new point is located using the history 
starting at the simplex of DT(Sk-i) containing the point. 

Analysis. In the analysis of the algorithm we will assume that the expected 
complexity of the Delaunay triangulation of a random sample of the given 
point set is linear in the size of the sample. This is the case for points in 
the plane, but it is also a realistic assumption for points sampled from a 
(d — l)-dimensional surface in IR d . The analysis can be adapted to the case 
where this assumption does not hold, yielding additional terms depending 
on the complexity of the triangulation. In the following theorem we bound 
the run-time in terms of the cost of computing the nearest-neighbor graph. 
Note that this bound holds for the standard real-RAM model and with no 
assumption about the spread of the point set. 

Theorem 1 Let P C IR d be a set of N points in general position such that 
the expected complexity of the Delaunay triangulation of a random sample R 
of P of size r is in 0(r). Algorithm^ constructs the Delaunay triangulation 
of P given in a random order in expected time F(N) + Y1T=\ ( m — i)F(2 l c) + 
O(N), where c is the (constant) size of the first round, m = |~log 2 (iV/c-|- 
and F{k) denotes the time needed to compute the nearest-neighbor graph of 
a subset of P of size k. 

Proof: We will analyse the cost of Step 3.3. Step 1 is only a conceptual 
step and Step 2 takes constant time. Step 3.1 takes constant time per loop (of 
Step 3). For a given k the nearest-neighbor graphs of T fc _iUi?fc-i, T k ^ 2 ^Rk-2, 
. . . , TiURi (or possibly fewer) are computed in Step 3.2. The size of these sets 
are bounded by 2 k c, 2 k ~ 1 c, . . . , 2c, where c = \Ri\ (except for k — 1 = m — 1, 
where \T k -i U Rk-i\ = N). Summing up over the loop of Step 3 this yields a 
cost of F(N) + Y,T=i ( m ~ i)F{2 i c) with m = \log 2 (N/c) + 1] . 

We now bound the cost of Step 3.3 for a given round (k > 1). The size of 
R k is c2 k ~ 1 . It suffices to prove that the cost of Step 3.3 is in 0(\R k \). For this 
we construct sets T' k _^ . . . ,Tq C P such that any point of R k has the same 
probability to be included into % U T[ (0 < % < k - 1). Let T' k _ x := 0. We 
construct T[ (k — 1 > % > 0) as follows: First we add each point of T[ +l with 
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probability 1/2. Second from each connected component in NNG(Tj +1 US , i+ i) 
with a vertex in Si+i we add each point of T i+ \ with probability 1/2. Note 
that the choices do not need to be independent. For all other connected 
components we add each point of T i+1 excluding the first and second (in the 
insertion order) with probability 1/2. Note that for a connected component 
all points have the same probability of being the first point of the component 
(in this case it is included in T; C TjUT/). Likewise all points have the same 
probability of being the second point (in this case it is not included), and 
likewise the same probability that it is one of the remaining points (in this 
case it is included with probability 1/2). Overall we get for all i that any 
point of R k is included into Tj U T[ with probability 2 l ~ k+1 and the expected 
size of Ti U T[ is |i4|2'~ fe+1 . 

We first bound the cost of locating a point p G Tj in DT(S,- + i). Since in 
the previous step we located p in DT(Sj), we can locate p using the history 
starting at a conflicting simplex of p in DT(Sj). Since Sj + i is a random 
subset of P and the points of Sj+i were inserted in a random order, the 
expected cost of locating p would be 0(\og(\Sj + i\/\Sj\)) = 0(1) if p were 
a random point of R k [11]. This is not the case, but for a random point 
of Tj U Tj it would be the case. The cost of locating all points of Tj in 
DT(Sj+i) is bounded by the cost of locating all points of T 3 -UT- in DT(Sj +1 ) 
(knowing a conflict in DT(Sj) for each point). The expected cost of this is 
in 0(E(Tj U Tj)) = 0(2i- k+1 \R k \). 

This gives us the expected cost of locating one conflicting simplex for each 
point p G Tj. We actually need to find the simplex in DT(Sj+i) containing 
the point p. This can be done by locally searching all conflicting simplices, 
one of which contains the point. The cost of this is therefore proportional to 
the number of conflicts a point has with DT(Sj +1 ). The cost for all points 
in Tj can again be bounded by the total number of conflicts of Tj U Tj with 
DT(S j+1 ) which is expected to be in 2^ k+1 0{\R k \). 

Second we bound the cost for locating the points of T J+ i in DT(5y + i). For 
locating a point of T J+1 , we traverse the connected components of NNG(T J+1 U 
Sj + \). For each component we start at a point for which we know the location 
in DT(Sj + i), i.e., a point from Sj+\UTj. Assume we traverse the edge between 
p and q where p is already located and q needs to be located. The point q is 
located by walking along the line segment pq, i.e., by traversing the Delaunay 
triangulation along pq. The cost corresponds to the number of intersected 
simplices. Any simplex intersected is either in conflict with p or with q, i.e., 
p or q lie in its circumsphere. If p G Sj + i we additionally have the cost of 
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searching for the <i-simplices adjacent to p that is the first simplex on the 
walk. 

Therefore the cost of the walk from p to q is bounded by the number of 
conflicts of q with simplices of DT(5^ + i) and - depending on whether p is 
in Tj +1 or Sj + i - by the number of conflicts of p or the number of faces at 
p. Any point can occur at most once as end point of a walk. Furthermore 
since the degree of NNG is in any fixed dimension bounded by the kissing 
number [10J, any point occurs only a constant times as starting point of a 
walk. Thus the total cost of walking is up to a constant factor bounded by 
the complexity of DT(Sj +1 ) and the expected total number of conflicts of 
Tj +1 with DT(S , j + i). By assumption the expected complexity of DT(S , J+ i) is 
linear in \Sj + i\ = (2 jf+1 — l)c. The expected number of conflicts of T J+ i with 
DT(Sj+i) we again bound by the expected number of conflicts of Tj +1 with 
DT(5 3+ i), which is 0(2 j ~ k+2 \R k \). Therefore the total expected cost of the 
round R k is in 0{^ k j l\2^ k+1 \R k \) = 0(\R k \). Summing up over all rounds 
yields an expected linear cost. ■ 

We now use that the nearest-neighbor graph of a point set with bounded 
spread can be computed in linear time. Note that the condition on the 
complexity of the Delaunay triangulation always holds in the plane, thus the 
algorithm computes the Delaunay triangulation of points in the plane with 
bounded spread in linear expected time. 
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Corollary 2 Let P C R be a set of N points in general position with 
bounded spread such that the expected complexity of the Delaunay triangula- 
tion of a random sample R C P of size r is in 0(r). Algorithm^ constructs 
the Delaunay triangulation of P given in a random order in linear expected 
time on a real-RAM with a constant-time fl,oor function restricted to logiV 
bits. 

We would like to note that the analysis can be extended to the case where 
the Delaunay hierarchy [12] is used instead of a history (the hierarchy is 
then built level by level) and also to the case of biased randomized insertion 
orders [TJ 0]. 
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